 /*
	Create 3 graphs, each showing statistics from a 
	range of starting years, through to end of PWT in 2019.
	First is average growth, second is volatility (SD), and
	third is persistence (growth corr in adjacent years).
*/


	use "$output/combined_data.dta", clear
	
*** Define income groups
	local data pwt
	summ `data' if ccode == "USA" & year == 1960 
	local usa = `r(mean)'
	summ pwt if year == 1960, det 
	local lowthresh = `r(p25)'/`usa'
	local highthresh = `r(p75)'/`usa'
	gen us = `data' if ccode == "USA"
	sort year us
	by year: replace us = us[_n-1] if missing(us) & !missing(us[_n-1])
	gen catl = `lowthresh'*us
	gen catu = `highthresh'*us
	gen cat = "H" if `data' >= catu & !missing(`data')
	replace cat = "M" if `data' < catu & `data' >=catl & !missing(`data')
	replace cat = "L" if `data' < catl 
			
*** Time set
	egen ccodenum = group(ccode)
	tsset ccodenum year
	gen logy = log(pwt)
	gen g1 = (log(pwt/l1.pwt))*100
	gen sd = .
	gen corr = .

*** Step 1. Create country-level statistics over time

*** Growth
	gen finalyear_ = pwt if year == 2019
	bys ccode: egen finalyear = max(finalyear_)
	drop finalyear_
	gen growth = (log(finalyear/pwt)/(2019- year))*100

*** Volatility and persistence
	levelsof ccode, local(countries)
	forval year = 1960/2010 {
		foreach country of local countries {
			qui count if !missing(g1) & ccode == "`country'" & year == `year'
			if `r(N)' == 1 {
				summ g1 if ccode == "`country'" & year >= `year'
				replace sd = `r(sd)' if year == `year' & ccode == "`country'"	
				sort ccodenum year
				gen temp = l1.g1
				corr g1 temp if ccode == "`country'" & year >= `year'
				replace corr = `r(rho)' if year == `year' & ccode == "`country'"
				drop temp 
			}
		}
	}
		

*** Step 2. Aggregate those results to country income-groups 

	foreach type in H M L{
		gen mean_`type' = .
		gen sd_`type' = .
		gen lb_`type' = .
		gen ub_`type' = .
		gen corr_`type' = .
		gen corrlb_`type' = .
		gen corrub_`type' = .
		gen p95_`type' = .
		gen p5_`type' = .
	}

	forval year = 1960/2010 {
		foreach type in H M L {	
		
		*** Averages
			*summ growth if cat == "`type'" & year >= `year', det
			ci means growth if cat == "`type'" & year == `year'
			replace mean_`type' = `r(mean)' if year == `year'
			*replace p95_`type'= `r(p95)' if year == `year'
			*replace p5_`type'= `r(p5)' if year == `year'
			replace p95_`type'= `r(lb)' if year == `year'
			replace p5_`type'= `r(ub)' if year == `year'

		*** Volatility
			ci means sd if cat == "`type'" & year == `year'
			replace sd_`type' = `r(mean)' if cat == "`type'" & year == `year'
			replace lb_`type' = `r(lb)' if cat == "`type'" & year == `year'
			replace ub_`type' = `r(ub)' if cat == "`type'" & year == `year'
			
		*** Persistence
			ci means corr if cat == "`type'" & year == `year'
			replace corr_`type' = `r(mean)' if cat == "`type'" & year == `year'
			replace corrlb_`type' = `r(lb)' if cat == "`type'" & year == `year'
			replace corrub_`type' = `r(ub)' if cat == "`type'" & year == `year'
		}
	}
	
	gen year2 =year + .2
	gen year3 = year + .4
	
	save "$output/for_triptych.dta", replace

*** Step 3. Make graphs

	use "$output/for_triptych.dta", clear

	global hcol plb1 
	global mcol black
	global lcol plg1 
	
	keep year* *_H *_M *_L 
	duplicates drop
	foreach var of varlist *_H *_L *_M {
		sort year `var'
		by year: replace `var' = `var'[_n-1] if missing(`var') & !missing(`var'[_n-1])
	}
	duplicates drop
	#delimit ;
	
	*** Averages ;
		tw 	(rcap p95_H p5_H year if inrange(year, 1960, 2010),  lcolor(${hcol}%30)) 
			(rcap p95_M p5_M year2 if inrange(year, 1960, 2010), lcolor(${mcol}%30)) 
			(rcap p95_L p5_L year3 if inrange(year, 1960, 2010), lcolor(${lcol}%30)) 
			(sc mean_H year if inrange(year, 1960, 2010),  mcolor(${hcol}%70)) 
			(sc mean_M year2 if inrange(year, 1960, 2010), mcolor(${mcol}%70)) 
			(sc mean_L year3 if inrange(year, 1960, 2010), mcolor(${lcol}%70))
			,
			legend(cols(1) order(4 "High income" 5 "Middle income" 6 "Low income") region(color(none)) bmargin(small) ring(1) pos(3)) 
			plotregion(style(none) lcolor(none))  
			xlabel(1960(5)2010, angle(0)) 
			ylabel(1(1)5, glcolor(gs10))
			graphregion(fcol(white) lcol(white)) 
			ytitle("{&mu}", orientation(horizontal) size(large)) 
			xtitle("Initial Year") 
			xsize(8) ysize(3) scale(2)
		;

		graph export "$figures/growthdist_hml_sc.pdf", replace 
		;
		
	*** Persistence ; 
		tw 	(rcap corrlb_H corrub_H year  if inrange(year, 1960, 2010), lcolor(${hcol}%30)) 
			(rcap corrlb_M corrub_M year2 if inrange(year, 1960, 2010), lcolor(${mcol}%30)) 
			(rcap corrlb_L corrub_L year3 if inrange(year, 1960, 2010), lcolor(${lcol}%30)) 
			(sc corr_H year  if inrange(year, 1960, 2010), mcolor(${hcol}%70)) 
			(sc corr_M year2 if inrange(year, 1960, 2010), mcolor(${mcol}%70)) 
			(sc corr_L year3 if inrange(year, 1960, 2010), mcolor(${lcol}%70)), 
			legend(cols(1) order(4 "High income" 5 "Middle income" 6 "Low income") region(color(none)) bmargin(small) ring(1) pos(3)) 
			plotregion(style(none) lcolor(none))  
			xlabel(1960(5)2010, angle(0)) 
			ylabel(-.1(.1).4, glcolor(gs10))
			graphregion(fcol(white) lcol(white)) 
			ytitle("{&rho}", orientation(horizontal) size(large)) 
			xtitle("Initial Year") 
			xsize(8) ysize(3) scale(2)
		;
		graph export "$figures/corr_anngrowth_hml_sc.pdf", replace 
		;
		
	*** Volatility ;
		tw 	(rcap lb_H ub_H year  if inrange(year, 1960, 2010), lcolor(${hcol}%30)) 
			(rcap lb_M ub_M year2 if inrange(year, 1960, 2010), lcolor(${mcol}%30)) 
			(rcap lb_L ub_L year3 if inrange(year, 1960, 2010), lcolor(${lcol}%30)) 
			(sc sd_H year  if inrange(year, 1960, 2010), mcolor(${hcol}%70)) 
			(sc sd_M year2 if inrange(year, 1960, 2010), mcolor(${mcol}%70)) 
			(sc sd_L year3 if inrange(year, 1960, 2010), mcolor(${lcol}%70)), 
			plotregion(style(none) lcolor(none))  
			xlabel(1960(5)2010, angle(0)) 
			ylabel(3(1)9, glcolor(gs10))
			graphregion(fcol(white) lcol(white)) 
			ytitle("{&sigma}", orientation(horizontal) size(large)) 
			xtitle("Initial Year")  
			legend(order(4 "High Income" 5 "Middle Income" 6 "Low Income") region(color(none)) ring(1) bmargin(small) cols(1) pos(3))
			xsize(8) ysize(3) scale(2)
		;
		graph export "$figures/stdev_anngrowth_hml_sc.pdf", replace 
		;
	#delimit cr
	


****************************
*** FIXED SAMPLE VERSION ***
****************************

	use "$output/combined_data.dta", clear
	
*** Define income groups
	local data pwt
	summ `data' if ccode == "USA" & year == 1960 
	local usa = `r(mean)'
	summ pwt if year == 1960, det 
	local lowthresh = `r(p25)'/`usa'
	local highthresh = `r(p75)'/`usa'
	gen us = `data' if ccode == "USA"
	sort year us
	by year: replace us = us[_n-1] if missing(us) & !missing(us[_n-1])
	gen catl = `lowthresh'*us
	gen catu = `highthresh'*us
	gen cat = "H" if `data' >= catu & !missing(`data')
	replace cat = "M" if `data' < catu & `data' >=catl & !missing(`data')
	replace cat = "L" if `data' < catl 
	
	keep if year >= 1980 
	bys ccode: egen count = count(pwt)
	keep if count == 40 
	drop count 

			
*** Time set
	egen ccodenum = group(ccode)
	tsset ccodenum year
	gen logy = log(pwt)
	gen g1 = (log(pwt/l1.pwt))*100
	gen sd = .
	gen corr = .

*** Step 1. Create country-level statistics over time

*** Growth
	gen finalyear_ = pwt if year == 2019
	bys ccode: egen finalyear = max(finalyear_)
	drop finalyear_
	gen growth = (log(finalyear/pwt)/(2019- year))*100

*** Volatility and persistence
	levelsof ccode, local(countries)
	forval year = 1980/2010 {
		foreach country of local countries {
			qui count if !missing(g1) & ccode == "`country'" & year == `year'
			if `r(N)' == 1 {
				summ g1 if ccode == "`country'" & year >= `year'
				replace sd = `r(sd)' if year == `year' & ccode == "`country'"	
				sort ccodenum year
				gen temp = l1.g1
				corr g1 temp if ccode == "`country'" & year >= `year'
				replace corr = `r(rho)' if year == `year' & ccode == "`country'"
				drop temp 
			}
		}
	}
		

*** Step 2. Aggregate those results to country income-groups 

	foreach type in H M L{
		gen mean_`type' = .
		gen sd_`type' = .
		gen lb_`type' = .
		gen ub_`type' = .
		gen corr_`type' = .
		gen corrlb_`type' = .
		gen corrub_`type' = .
		gen p95_`type' = .
		gen p5_`type' = .
	}

	forval year = 1980/2010 {
		foreach type in H M L {	
		
		*** Averages
			*summ growth if cat == "`type'" & year >= `year', det
			ci means growth if cat == "`type'" & year == `year'
			replace mean_`type' = `r(mean)' if year == `year'
			*replace p95_`type'= `r(p95)' if year == `year'
			*replace p5_`type'= `r(p5)' if year == `year'
			replace p95_`type'= `r(lb)' if year == `year'
			replace p5_`type'= `r(ub)' if year == `year'

		*** Volatility
			ci means sd if cat == "`type'" & year == `year'
			replace sd_`type' = `r(mean)' if cat == "`type'" & year == `year'
			replace lb_`type' = `r(lb)' if cat == "`type'" & year == `year'
			replace ub_`type' = `r(ub)' if cat == "`type'" & year == `year'
			
		*** Persistence
			ci means corr if cat == "`type'" & year == `year'
			replace corr_`type' = `r(mean)' if cat == "`type'" & year == `year'
			replace corrlb_`type' = `r(lb)' if cat == "`type'" & year == `year'
			replace corrub_`type' = `r(ub)' if cat == "`type'" & year == `year'
		}
	}
	
	gen year2 =year + .2
	gen year3 = year + .4
	
	save "$output/for_triptych_fixed.dta", replace

*** Step 3. Make graphs

	use "$output/for_triptych_fixed.dta", clear

	global hcol plb1 
	global mcol black
	global lcol plg1 
	
	keep year* *_H *_M *_L 
	duplicates drop
	foreach var of varlist *_H *_L *_M {
		sort year `var'
		by year: replace `var' = `var'[_n-1] if missing(`var') & !missing(`var'[_n-1])
	}
	duplicates drop
	#delimit ;
	
	*** Averages ;
		tw 	(rcap p95_H p5_H year if inrange(year, 1980, 2010),  lcolor(${hcol}%30)) 
			(rcap p95_M p5_M year2 if inrange(year, 1980, 2010), lcolor(${mcol}%30)) 
			(rcap p95_L p5_L year3 if inrange(year, 1980, 2010), lcolor(${lcol}%30)) 
			(sc mean_H year if inrange(year, 1980, 2010),  mcolor(${hcol}%70)) 
			(sc mean_M year2 if inrange(year, 1980, 2010), mcolor(${mcol}%70)) 
			(sc mean_L year3 if inrange(year, 1980, 2010), mcolor(${lcol}%70))
			,
			legend(cols(1) order(4 "High income" 5 "Middle income" 6 "Low income") region(color(none)) bmargin(small) ring(1) pos(3)) 
			plotregion(style(none) lcolor(none))  
			xlabel(1980(5)2010, angle(0)) 
			ylabel(1(1)5, glcolor(gs10))
			graphregion(fcol(white) lcol(white)) 
			ytitle("{&mu}", orientation(horizontal) size(large)) 
			xtitle("Initial Year") 
			xsize(8) ysize(3) scale(2)
		;

		graph export "$figures/growthdist_hml_sc_fixed.pdf", replace 
		;
		
	*** Persistence ; 
		tw 	(rcap corrlb_H corrub_H year  if inrange(year, 1980, 2010), lcolor(${hcol}%30)) 
			(rcap corrlb_M corrub_M year2 if inrange(year, 1980, 2010), lcolor(${mcol}%30)) 
			(rcap corrlb_L corrub_L year3 if inrange(year, 1980, 2010), lcolor(${lcol}%30)) 
			(sc corr_H year  if inrange(year, 1980, 2010), mcolor(${hcol}%70)) 
			(sc corr_M year2 if inrange(year, 1980, 2010), mcolor(${mcol}%70)) 
			(sc corr_L year3 if inrange(year, 1980, 2010), mcolor(${lcol}%70)), 
			legend(cols(1) order(4 "High income" 5 "Middle income" 6 "Low income") region(color(none)) bmargin(small) ring(1) pos(3)) 
			plotregion(style(none) lcolor(none))  
			xlabel(1980(5)2010, angle(0)) 
			ylabel(-.1(.1).4, glcolor(gs10))
			graphregion(fcol(white) lcol(white)) 
			ytitle("{&rho}", orientation(horizontal) size(large)) 
			xtitle("Initial Year") 
			xsize(8) ysize(3) scale(2)
		;
		graph export "$figures/corr_anngrowth_hml_sc_fixed.pdf", replace 
		;
		
	*** Volatility ;
		tw 	(rcap lb_H ub_H year  if inrange(year, 1980, 2010), lcolor(${hcol}%30)) 
			(rcap lb_M ub_M year2 if inrange(year, 1980, 2010), lcolor(${mcol}%30)) 
			(rcap lb_L ub_L year3 if inrange(year, 1980, 2010), lcolor(${lcol}%30)) 
			(sc sd_H year  if inrange(year, 1980, 2010), mcolor(${hcol}%70)) 
			(sc sd_M year2 if inrange(year, 1980, 2010), mcolor(${mcol}%70)) 
			(sc sd_L year3 if inrange(year, 1980, 2010), mcolor(${lcol}%70)), 
			plotregion(style(none) lcolor(none))  
			xlabel(1980(5)2010, angle(0)) 
			ylabel(3(1)9, glcolor(gs10))
			graphregion(fcol(white) lcol(white)) 
			ytitle("{&sigma}", orientation(horizontal) size(large)) 
			xtitle("Initial Year")  
			legend(order(4 "High Income" 5 "Middle Income" 6 "Low Income") region(color(none)) ring(1) bmargin(small) cols(1) pos(3))
			xsize(8) ysize(3) scale(2)
		;
		graph export "$figures/stdev_anngrowth_hml_sc_fixed.pdf", replace 
		;
	#delimit cr
	
exit
